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ABSTRACT 


The  thermal  fission  of  Uranium-235  has  been  studied  by  developin; 
a  computer  code  that  uses  evaporation  theory  based  on  the  Weisskopf 
model  for  neutron  emission  from  individual  fission  fragments.  This 
computer  code  calculates  a  prompt  neutron  spectrum  for  the  thermal 
fission  of  Uranium-235,  calculates  the  total  gamma  ray  energy  for  the 
fission  of  Uranium-235,  and  calculates  the  average  prompt  neutron  kin¬ 
etic  energy.  Until  now,  evaporation  theory  codes  used  average  values 
of  mass  and  energy  in  determining  prompt  neutron  spectra,  and  this 
code  was  developed  to  use  all  possible  fission  fragment  combinations 
In  order  to  eliminate  this  averaging.  This  approach  Is  necessary  to 
compute  certain  measured  data  such  as  the  probability  that  0,1,2,... 

'Ui.  C.S 

neutrons  are  emitted  per  fission.  This  codeAutil1zes  subroutines  to 
calculate  rest  masses,  corrects  the  calculated  rest  mass  to  agree  with 
experiment,  and  to  determine  the  emitted  neutron  kinetic  energies. 

Evaporation  theory  codes  have  predicted  low  total  gamma  energy 
release  In  fission,  but  agree  quite  favorably  with  experiments  in  pre¬ 
dicting  the  averge  prompt  neutron  kinetic  energy.  The  prompt  neutron 
spectrum  obtained  in  this  code  has  a  high  tall  region,  that  is,  below 
0.5  MeV.  The  average  kinetic  energy  of  the  neutrons  agree  very  well 
with  experiment,  and  the  total  gamma  energy  Is  lower  than  experiment, 
as  previous  evaporation  theory  codes  have  been.  f 


CHAPTER  1 


INTRODUCTION  AND  HISTORICAL  PERSPECTIVE 

Fission  is  the  process  which  occurs  when  a  fissionable  nucleus 
captures  a  neutron  and  splits  into  two  lighter  nuclei,  increasing  the 
stability  of  the  system.  When  fission  occurs,  up  to  a  half-dozen 
neutrons  are  emitted  along  with  gamma  rays,  and  at  some  time  later, 
beta  particles  and  neutrinos.  The  fission  fragments  forme'  as  a  re¬ 
sult  of  fission  are  neutron  rich  and  must  decay  to  a  more  stable  con¬ 
figuration.  One  of  the  most  Important  products  of  the  fission  proc¬ 
ess  is  the  prompt  neutrons  which  are  emitted  within  the  first  10“17 
seconds  of  the  split  of  the  parent  nucleus. 

In  the  design  of  nuclear  reactors,  it  Is  Important  to  know  the 
energy  spectrum  for  the  neutrons  emitted  as  a  result  of  the  fission¬ 
ing  of  the  compound  nucleus  and  the  energy  of  the  gamma  rays  emitted 
from  the  fission  process.  The  neutrons  are  necessary  for  the  sus¬ 
tainment  of  the  fission  process,  and  the  gamma  rays  are  a  biological 
hazard  and  require  shielding. 

James  Terrell,  in  1958,  published  his  study  entitled  "Fission 
Neutron  Spectra  and  Nuclear  Temperatures".  His  paper  considered  "the 
agreement  with  experiment  of  various  theoretical  predictions  as  to 
the  energy  spectrum  of  fission  neutrons".1  He  also  gathered  Informa¬ 
tion  concerning  nuclear  temperatures  from  the  measured  spectra,  and 
looked  at  how  the  neutron  multiplicity  (that  Is,  the  average  number 
of  neutrons  emitted  in  the  fissioning  process)  varies  with  the  exci¬ 
tation  energy.  He  also  calculated  a  center  of  mass  neutron  spectrum 


and  calculated  the  effect  of  anisotropy  of  neutron  emission  on  the 
laboratory  spectrum.  Terrell  concludes  that  all  experimental  energy 
distributions  for  fission  neutrons  are  indistinguishable  from 
Maxwelliam  distributions  and  shows  that  distributions  of  essentially 
this  form  are  predicted  by  Weisskopf's  evaporation  theory. 
"Calculation  from  evaporation  theory  gives  the  observed  distribution 
of  neutron  numbers;  but  calculated  total  gamma-ray  energies  are  too 
low,  unless  gamma  rays  are  emitted  much  more  rapidly  than  commonly 
assumed."!  Neutron  multiplicity  in  Terrell's  paper  is  determined  to 
be  2.46  neutrons/fission. 

In  Terrell's  paper,  calculations  are  performed  using  an  average 
light  mass  and  an  average  heavy  mass.  His  average  mass  of  the  heavy 
fragment  to  the  light  fragment  is  1.47  for  Uranium  235.  A  compila¬ 
tion  of  the  average  energy  of  fission  neutrons  is  presented  in  his 
paper  based  on  experimental  methods.  They  are  for  U-235  +  neutron: 


Experimenter! 

Energy  Range  (MeV) 

T  (MeV) 
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0.18  -  3 
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Time  of 
Flight 

Frye  4  Rosen 

0.35  -  12 

1.9806 

Photoplate 

Mulkhin 

1  -  12 

— 

Photoplate 

Nereson 

0.4  -  7 

— 

Photoplate 

Watt 

i 

3.3  -  17 

2.00 

Proportional 

Counter 

The  more  recent  research  on  the  prompt  neutron  energy  spectrum 
has  been  performed  by  David  G.  Madland  and  J.  Rayford  Nix,  at  Los 
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Alamos  Scientific  Laboratory.  Agnin,  average  values  of  the  light  and 
heavy  fission  fragments  are  used.  They  use  standard  evaporation 
theory  to  calculate  the  neutron  energy  spectrum  in  the  center-of-mass 
of  a  given  fission  fragment,  and  then  transform  into  the  laboratory 
system,  taking  into  account  that  the  average  velocity  of  the  light 
fragment  is  higher  than  the  average  velocity  of  the  heavy  fragment. 
The  energy  released  in  fission  is  based  on  seven  heavy  fission  frag¬ 
ments.  These  fragments  are  2  Iodine  Isotopes,  3  Xenon  Isotopes  and  2 
Cesium  isotopes.  For  the  total  average  fission-fragment  kinetic 
energies,  Madland  and  Nix  used  the  experimental  results  of  Un'.k  et 
al.,  which  are  appropriate  to  thermal  neutron  induced  fission,  how¬ 
ever,  these  are  average  values  only.  Their  results  show  an  average 

neutron  multiplicity  of  2.403  for  thermal  neutron  Induced  fission  of 

* 

Uranium-235.  In  the  laboratory  system,  the  average  energy  of  the 
prompt  neutrons  for  0.53  MeV  Induced  fission  and  0.60  MeV  Induced 
fission  Is  2.046  and  2.048  respectively.  No  kinetic  energy  results 
for  thermal  fission  are  given  In  their  paper. 

As  seen  In  the  literature,  there  Is  a  need  to  apply  first  prin¬ 
ciples  to  the  fission  process  for  all  possible  fragment  combinations 
In  order  to  remove  the  averaging  of  fragment  weight  and  kinetic 
energy.  The  code  developed  In  this  thesis  uses  first  principles  on 
all  fission  fragment  combinations  for  the  thermally  induced  fission 
of  Uranium  235.  Major  Items  of  Interest  are  (1)  to  develop  a  neutron 
energy  spectrum  for  thermal  fission  of  Uranium-235,  (2)  determine  the 
neutron  multiplicity  of  Uranium-235,  (3)  determine  the  total  gamma 
energy  available  on  the  fragments.  Total  gamma  energy  Is  determined 


by  the  initial  excitation  energy  in  the  case  when  no  neutrons  are 
emitted,  and  the  residual  energy  in  the  case  when  neutrons  are  emit¬ 
ted,  but  energy  remains  which  is  below  that  required  to  emit  a  neu¬ 
tron,  This  total  gamma  energy  is  weighted  by  the  yield  of  the 
particular  nuclide  and  displayed  as  a  result  in  the  code. 

Another  development  which  this  code  will  allow  is  the  determina¬ 
tion  of  the  number  of  neutrons  emitted  from  each  possible  fragment. 
This  allows  one  to  compute,  by  knowing  the  fractional  yield  of  the 
nuclide  and  the  number  of  neutrons  emitted  by  the  fragment,  the  con¬ 
tribution  of  the  neutron  multiplicity  from  each  individual  fragment. 

The  description  of  fission  theory  in  Chapter  2  is  presented 
largely  in  the  form  of  quotations.  The  application  to  the  code  de¬ 
velopment  is  presented  in  parallel  to  the  theory.  In  brief,  the 
yield  of  each  primary  fission  fragment  is  determined  as  the  measured 
A-chain  yield  multiplied  by  a  Gaussian  probability  distribution  in 
charge.  For  a  given  primary  fission  fragment  pair  the  rest  mass 
energies  are  computed  and  subtracted  from  the  mass  energy  of  tne 
fissioning  nucleus.  This  energy  decrement  appears  as  the  sum  of  the 
fragment  excitation  energies  plus  the  kinetic  energies  of  the  sepa¬ 
rating  fragments.  The  latter  Is  determined  in  ten^  of  the  fragment 
charges  and  the  distance  RR  between  the  fragments  at  the  instant  of 
scission.  For  an  assumed  RR  the  fragment  kinetic  energies  are  aeter- 
mined  and,  from  the  sum  of  excitation  energies,  a  plausible  assump¬ 
tion  results  in  the  initial  excitation  energy  on  each  primary  fission 
fragment.  Next,  computing  the  mass  of  the  nucleus  which  would  result 
if  neutron  emission  occurred,  it  is  determined  whether  neutron 
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emission  is  energetically  possible.  If  so,  the  neutron  kinetic 
energy  in  the  fission  fragment  frame  of  reference  is  determined  by 
evaporation  theory,  and  we  are  again  left  with  an  excited  nucleus. 

The  process  is  repeated  until  neutron  emission  is  not  energetically 
possible,  after  which  the  total  gamma  deexcitation  energy  is  record¬ 
ed;  the  gamma  emission  spectrum  is  not  computed.  Finally,  the  neu¬ 
tron  kinetic  energy  is  transformed  to  the  laboratory  coordinate 
system.  The  calculations  are  repeated  for  all  fission  fragment  pairs 
and  for  various  values  of  RR.  The  best  value  of  RR  is  chosen  to  give 
the  best  agreement  with  measured  data. 


CHAPTER  2 


THEORY 


2.1  Mass  Distribution 

"One  of  the  most  characteristic  features  of  fission  is  the 
asymmetric  division  of  the  fissioning  nucleus,  and  for  many  years  our 
most  complete  knowledge  of  the  mass  division  came  from  radiochemical 
research. "2  "Chat  particular  fragment  nuclides  are  produced  by  the 
given  nucleus  is  a  matter  of  chance  and  the  range  of  gross  fission 
products  Is  roughly  from  bromine  to  barium  in  the  periodic  table. 

The  concentration  of  fission  nuclides  depends  on  the  atomic  mass  and 
the  distribution  curve  has  a  curious  saddleback  shape  (see  Fig.  2.1). 
These  are  two  well-defined  maxima,  at  A  a  95  and  A  *  140  (for  the 
fission  of  U-235) ,  roughly.  The  total  yield  Is  200%  since  there  are 
two  fragments  per  fission.  The  total  number  of  Identified  fission 
nuclides  Is  about  300,  Including  nearly  200  different  0"  emitters."-* 

"The  first  work  on  fission  yield  and  the  Introduction  of  the 
concept  of  fission  yield  were  due  to  Fermi  and  his  co-workers  at 
Columbia."2 

"The  determination  of  the  fission  yield  of  a  specific  species 
consists  of  a  number  of  steps. 

1.  A  measured  amount  of  non-radloactive  carrier  material  of  a 
given  fission  product  element  is  added  to  a  solution  of 
uranium  In  which  a  known  number  of  fission  events  has 


occurred 


2 


.  If  it  is  necessary,  chemical  treatment  is  given  this 

solution  to  insure  complete  isotopic  exchange  of  the  stable 
and  radioactive  isotopes  of  the  element. 

3.  The  solution  is  subjected  to  an  analytical  procedure  to 
separate  the  element  from  the  solution  in  a  state  of 
chemical  and  radioactive  purity. 

4.  The  fractional  recovery  of  the  inert  carrier  is  determined 
by  some  quantitative  analytical  method.  The  chemical 
recovery  of  the  tracer  element  is  assumed  equal  to  that  for 
the  inert  carrier  material. 

5.  The  radiation  of  the  purified  radioelements  are  measured  to 
identify  the  isotopes  and  to  determine  the  absolute  amounts 
of  each  species.  Correction  is  made  for  radioactive  decay 
from  the  time  of  fission  to  the  time  of  counting. 

6.  From  the  counting  data,  the  chemical  yield  data,  and  the 
known  number  of  fission  events  the  fission  yield  Is 
calculated.  The  fission  yield  is  defined  as  the  percentage 
of  fissions  leading  to  the  formation  of  a  measured  product. 

It  Is  to  be  noted  that  the  radiochemical  results  do  not  In 
general  give  the  independent  yield  of  the  specific  isotope  measured. 
Usually,  the  experimentally  determined  yield  is  the  cumulative  yield 
of  the  specific  Isotope,  including  any  precursors  which  have 
undergone  decay  to  the  specific  isotope  before  the  chemical  Isolation 
occur red. "2 

"The  largest  error  in  determination  of  a  fission  yield  lies  in 
the  evaluation  of  the  absolute  disintegration  rate  of  the  radioactive 


species.  In  order  to  correct  for  the  effects  of  absorbers  [such  as 
the  window  of  the  detector,  the  material,  the  source  covering,  etc.], 
it  is  usual  to  establish  empirically  the  initial  half-thickness  of 
the  radiation  emitted  by  the  source  by  means  of  the  determination  of 
an  absorption  curve.  In  many  instances  such  a  curve  is  complex.  The 
establishment  of  an  initial  half-thickness,  then,  becomes  a  matter  of 
interpretation  and  individual  judgment".4 

"In  thermal  neutron  fission  of  most  nuclei,  an  asymmetric  mass 
split,  with  a  ratio  of  heavy-to-light  fragment  mass  (M^/M(_)  of  about 
1.4,  is  much  more  likely  than  a  symmetric  mass  split.  With 
increasing  bombarding  energy  the  valley  between  the  two  humps  in  the 
fission-yield  curve  gradually  fills  in,  so  that,  for  example,  with  14 
MeV  neutrons  on  U-235,  the  peak  to  valley  ratio  is  only  about  6 
(compared  to  about  600  for  thermal  neutrons).  When  the  bombarding 
energy  reaches  about  50  meV,  the  fission-yield  curves  for  the  highly 
fissile  elements  (Z  >  90)  exhibit  single  bread  humps  with  the  valley 
completely  filled  In.  It  has  recently  been  established  that  for  the 
very  heaviest  elements,  that  is  in  the  region  of  fermium,  the 
probability  of  a  symmetrical  mass  split  again  Increases  relative  to 
the  asymmetric  modes,  until,  for  the  thermal -neutron  fission  of 
Fm-257,  symmetric  fission  predoml nates. "5 

2.2  Charge  Distribution 

"An  important  part  of  the  information  that  one  would  like  to 
have  about  the  fission  process  is  the  division  of  nuclear  charge 
between  the  primary  fission  fragments.  Unfortunately,  to  determine 
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this  is  a  difficult  experimental  problem  and  the  available  data  are 
limited.  The  reason  for  this  difficulty  is  that  the  primary 
fragments  are  so  far  from  beta  stability  that  most  of  them  have  very 
short  half-lives.  Hence,  by  the  time  the  necessary  chemical 
separations  have  been  carried  out,  the  primary  products  have  been 
completely  converted  into  different  elements.  This  is  not  true  in 
the  case  of  shielded  nuclides,  and  their  fission  yields  are  of 
necessity  independent  rather  than  cumulative  chain  yields.  A 
shielded  nuclide  is  one  which  cannot  be  formed  by  beta  decay  because 
the  isobaric  nuclide  of  the  next  lower  atomic  number  is  stable. "2 
"Several  postulates  have  been  made  concerning  the  manner  in 
which  the  nuclear  charge  of  the  fissioning  nucleus  distributes 
itself  between  the  two  fragments: 

(1)  The  unchanged  charge  distribution  (UCD)  or  constant  charge 
ratio  (CCR)  postulate  predicts  that  the  two  fragments 
will  have  he  same  charge  density  (Z0/A0)  as  Lhe 
fissioning  nucleus.  This  type  of  charge  division  leads  to 
chain  lengths  which  are  longer  in  the  light  fragment. 
Approximately  equal  chain  lengths  are  observed  in 
low-energy  fission. 

(2)  The  postulate  of  minimum  potential  energy  (MPE)  proposed  by 
R.D.  Present  (1947),  postulated  that  nuclear  charge  dis¬ 
tributes  itself  between  fragments  such  that  a  minimum  is 
achieved  between  the  nuclear  potential  energy  and  the 
Coulomb  energy.  This  corresponds  to  maximizing  the 
excitation  energy.  The  postulate  predicts  chain  lengths 


of  the  light  fragments  which  are  shorter  than  those  of  the 
heavy  fragments. 

(3)  The  postulate  of  maximum  energy  release  (MER)  was  proposed 
by  Kennett  and  Thode  (1956)  to  explain  better  the  primary 
yields  of  nuclei  in  the  region  of  Z=50.  The  maximum 
energy  release  is  calculated  as  the  difference  between  the 
mass  of  the  fissioning  nucleus  and  the  masses  of  the 
primary  fission  products  and  neutrons  emitted.  The  most 
probable  nuclear  charge  which  results  from  this  method 
depends  upon  the  mass  formula  that  is  used  to  determine 
the  primary  fission  product  masses. 

(4)  The  postulate  of  equal  charge  displacement  (ECD)  is 
essentially  an  empirical  correlation  of  low-energy 
experimental  data.  The  two  basic  assumptions  of  this 
postulate  are  (a)  the  most  probable  charges  for 
complementary  fission  products  be  an  equal  number  of 
charge  units  away  from  the  valley  of  beta  stability  and 
(b)  the  distribution  function  about  the  most  probable 
charge  is  symmetric  and  the  same  for  all  mass  splits."6 

"The  postulate  that  the  most  probable  division  of  charge  in 
fission  is  the  division  that  gives  rise  to  fragments  of  equal 
effective  chain  length  (equal  charge  displacement)  may  now  be  taken 
as  a  working  hypothesis  of  the  variation  In  nuclear  charge  among  the 
primary  fission  products.  It  should  be  emphasized  that  this  is 
essentially  an  empirical  hypothesis.  Its  validity  may  be  based  on 
the  fact  that  there  is  competition  between  the  effects  of  splitting 


of  the  nucleus  to  give  a  charge  distribution  not  to  greatly  different 
from  that  of  the  fissioning  nucleus,  with  perhaps  slightly  asymetric 
charge  distribution  as  suggested  by  Present,  (Phys.  Rev,  72:7  (1947)) 
and  the  effect  of  splitting  to  give  minimal  potential  energy. "7 

"The  shape  of  the  distribution  function  which  represents  the 
relative  probability  P(Z)  of  forming  a  product  with  atomic  number  Z 
is  assumed  to  be  a  Gaussian  given  by  the  relationship: 

P(Z)  ■  [l/(CTr)i/2]  exp[-(Z-Zp)2/c]  (2.1) 

where  c  is  a  constant  empirically  determined."® 

"This  relationship  holds  for  all  mass  numbers,  with  the  param¬ 
eter  c  =  0.79±0.14  (for  23Su  fission  by  thermal  neutrons).  This 
corresponds  to  a  full  width  at  half  maximum  of  1.50+0.12  charge 
units.  (Wahl  et  al .  "Nuclear-Charge  Distribution  in  Low-Energy 
Fission,"  Physical  Review  126,  1112  (1962)). "5 

"From  the  assumption  of  equal  charge  displacement,  one  may  write 

<ZAL  '  V  ’  <ZAH  -  W  ‘Z-2> 

where  Z a  refers  to  the  most  stable  charge,  Zp  to  the  most  probably 
charge,  and  the  subscripts  L  and  H  refer  to  the  light  and  heavy 
complementary  primary  products,  respectively.  The  following 
relationships  are  valid  for  binary  fission: 

ZPL  +  ZPH  *  Z0  (2,3) 

and 

\  *  Ah  ■  A0  -  v  (2.4) 

where  v  is  the  number  of  neutrons  emitted  in  the  mass  split.  The 
most  probable  charge  of  a  fission  product  of  mass  A  is  given  by: 


Thus  the  Zp  function  is  defined  in  terms  of  the  most  stable 
masses."6 

Wahl  was  able  to  determine  the  primary  and  cumulative  yields  of 
nine  short  lived  isotopes  of  krypton  and  xenon  by  measuring  the 
fractional  cumulative  yields  of  inert  gases  produced  as  fission 
products  and  captured  in  a  thin  layer  of  barium  stearate  powder. 

From  this  work,  Wahl  developed  an  empirical  curve  for  Zp.  "Several 
of  Wahl’s  students  have  contributed  newer  data  on  independent  and 
cumulative  yields  which  are  significant  for  an  analysis  of  the  equal 
charge  displacement  hypothesis  in  its  various  formulations;  the 
comprehensive  summary  paper  on  this  work  is  a  key  paper  in  the 
literature  on  charge  distribution."2 

In  order  to  code  the  equal  charge  displacement  hypothesis  a  few 
variables  must  be  defined. 

AN(IA)  is  the  mass  of  the  IAth  nuclide. 

ZN(IA,IZ)  is  the  charge  number  of  the  IZtfl  nuclide  with  mass 
number  AN(IA). 

ZA(IA)  is  the  starting  Z  for  the  IAth  chain. 


YA(IA)/200.0  is  the  chain  yield  divided  by  200%,  since  the  yield 


is  given  by  thi  number  out  of  200  events  that  it  will  occur.  The 
independent  yield  is  given  as  the  product  of  the  chain  yield/200.0 
and  the  Probability  Distribution  for  the  specified  Z-chain. 

YZ( IA,IZ)  =  (YA(IA)/200.0)*C1*EXP(-((ZN(IA,IZ)-ZP(IA))**2)/CC) 

(2.7) 

After  all  yields  are  determined,  the  yield  weighted  average  mass 
may  be  determined  by  summing  up  all  of  the  products  of  yield  and  mass 


l  [( YA( IA)/200.00)*AN( IA)] 
IA 


(2.8) 


and  then  divided  by  the  sum  of  all  yields: 


YA( IA)/200.00 


(2.9) 


So,  the  yield  weighted  average  of  all  masses  can  be  written: 

l  (YA( IA)/200.00)*AN( I A) 

(M5T)  =  - fVA'MM  IVKR’hKS - 


(2.10) 


This  averaging  is  performed  for  two  groups.  First,  all  nuclides  with 
4  <  A  <  114  are  considered  In  the  light  group.  Second,  all  nuclides 
with  115  <  A  <  235  are  considered  In  the  heavy  group.  From  this 
averaging  process,  the  Average  heavy  fragment  and  average  light 
fragment  mass  may  be  determined. 


2.3  Hass  and  Binding  Energy 

"The  liquid-drop  model  provides  reasonable  explanations  for  many 
phenomena  which  are  inaccessible  to  the  shell  model.  These  phenomena 
are  mainly: 


(1)  Substantially  constant  density  of  nuclei  with  radius 
R  «  Ai/3. 

(2)  Systematic  dependence  of  the  neutron  excess  (N-Z)  on  A5/3 
for  stable  nuclides. 

(3)  Approximate  constancy  of  the  binding  energy  per  nucleon 
B/A,  as  well  as  its  small  but  definite  systematic  trends 
with  A  (see  Fig.  2.2) . 

(4)  Mass  differences  in  families  of  isobars  and  the  energies  of 
cascade  0  transitions. 

(5)  Systematic  variation  of  a  decay  energies  with  N  and  Z. 

(6)  Fission  by  thermal  neutrons  of  U-235  and  other  odd-N 
nuclides. 

(7)  Finite  upper  bound  Z  and  N  of  heavy  nuclides  produced  in 
reactions  and  the  non-existence  In  nature  of  nuclides 
heavier  than  [Pu-244]. 

which  involve  the  masses  and  binding  energy  of  nuclear  ground  levels; 
the  energetics  of  e-decay,  a  decay,  and  nuclear  reactions;  and  the 
energetics  of  nuclear  fission."8 

•  "The  basic  premise  of  the  liquid  drop  model  of  complex  nuclei  Is 
that  the  number  of  nucleons  in  a  typical  nucleus  is  sufficiently 
large  that  the  Individuality  of  nucleons  may  be  disregarded.  While 
the  treatment  of  such  a  system  is  most  rigorously  given  by  using  the 
methods  of  quantum  statistical  mechanics,  it  is  possible  to  arrive  at 
three  components  of  the  expression  for  the  energy  of  a  complex 
nucleus  by  simple  classical  considerations.  These  components  are  the 
volume  binding  energy,  the  surface  energy,  and  the  coulomb  energy. "9 


constant.  (From  R.D.  Evans,  The  Atomic  Nucleus.  P.  299) 


In  addition  to  the  three  terms  which  may  be  given  by  classical 
considerations,  there  are  also  additional  terms  in  the  semi -empirical 
mass  formula.  These  terms  are  called  the  symmetry  terms,  which 
Meyers  and  Swiatecki  include  in  the  volume  term,  the  correction  term 
to  the  Coulomb  energy  term  and  also  the  pairing  energy  term.  These 
terms  are  discussed  individually: 

(1)  Volume  Energy  Term:  "The  first  and  dominant  term, 
proportional  to  A  and  thus  to  the  nuclear  volume,  expresses 
the  fact  that  the  binding  energy  is  in  first  approximation 
proportional  to  the  number  of  nucleons.  This  Is  a  direct 
consequence  of  the  short  range  and  saturation  character  of 
the  nuclear  forces.  The  saturation  is  almost,  though  not 
entirely,  complete  when  four  nucleons,  two  protons  and  two 
neutrons.  Interact,  as  is  indicated  by  the  large  observed 
binding  energies  of  4He,  12C,  160  and  so  on."- 

(2)  Surface  Energy  Term:  "  nucleons  at  the  surface  'feel 
less  interaction'  than  nucleons  deep  inside;  compared  to 
the  latter,  the  external  medium,  vacuum,  deprives  them  of 
more  than  half  their  partners.  The  volume  energy  term 
over-estimates  the  total  binding,  since  the  nucleus  is 
finite.  The  lost  binding  is  proportional  to  the  number  of 
nucleons  at  the  surface,  that  is  to  the  surface  area.  By 
assuming  the  nuclear  drop  to  be  a  sphere.  Its  surface  area 
is  given  by  4*r0*A2/3  and  the  correction  to  be  applied  can 
be  written  as  UB)S  =  -asA2/3  where  as  is  a  constant."^ 
With  increasing  nuclear  size,  the  surface-to-volume  ratio 


decreases,  and  therefore  this  term  becomes  relatively 
unimoortant."5 

Coulomb  Energy  Term:  "The  above  constitute  the  main 
feature  of  an  electrically  neutral  drop.  We  must  now 
consider  the  Coulomb  interaction  which  acts  between 
protons.  This  is  a  repulsive  action,  and  therefore  its 
contribution  decreases  the  binding  energy.  It  is 
responsible  for  the  slow  decrease  of  B/A  beyond  A  2  60,  as 
Z  Increases  (Fig.  2.3).  For  all  light  nuclei,  this 
contribution  Is  small  since  the  strong  Interaction  is 
stronger  than  the  electromagnetic  Interaction.  If  this 
contribution  becomes  appreciable  for  high  Z,  it  is  because 
the  Coulomb  Interaction  is  of  'infinite'  range,  and 
therefore  acts  over  the  whole  volume,  whereas  the  strong 
Interaction  acts  over  a  region  of  radius  close  to  1.5 
Fermi s."10  "The  electrostatic  energy  of  a  uniformly 
charged  sphere  of  charge  q  and  radius  R  is  3/5  q2/R  and 
since  q  *  Ze  and  R  *  r0A1/3  for  a  nucleus  of  radius  R  and 
atomic  number  Z,  we  can  write  its  electrostatic  energy  as 
(3e2/5r0)Z2A-i/3.  The  coefficient  (3e2/5r0)  *  0.717  MeV 
corresponds  to  an  r0  value  of  1.205  Fermi s.  Because  of  its 
Z2  dependence,  the  coulomb  energy  becomes  Increasingly 
Important  as  Z  Increases  and  accounts  for  the  fact  that  all 
stable  nuclei  with  Z>20  contain  more  neutrons  than  protons 
(see  Fig.  3.4)  despite  the  symmetry  energy  that  maximizes 
nuclear  binding  for  N»Z."5 
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Fig.  2.3  Summary  of  the  semi -empirical  liquid-drop-model 
treatment  of  the  average  binding  energy  curve. 
Note  how  the  decrease  In  surface  energy  and  the 
Increase  In  coulomb  energy  cons  1  pi  re  to  produce 
the  maximum  observed  In  B/A  at  A~60.  (From 
R.O.  Evans,  The  Atomic  Nucleus,  p.  382) 


charged,  but  have  charge  distributions  with  diffuse 
boundaries.  The  diffuse  boundary  gives  rise  to  a 
correction  to  coulomb  energy  (lowering  it),  and  this  is 
expressed  in  the  fourth  term  of  the  semiempirical  mass 
formula. "5  This  correction  term,  introduced  by  Meyers  and 
Swiatecki  is  written  as:  (ir2e2/(2r0))  •  (d/r0)2Z2/A,  where 
d  is  given  by  0.5461  Fermis,  based  on  the  Stanford  electron 
scattering  experiments.  So  this  correction  term  can  be 
written:  1.211  Z2/A. 

(5)  Symmetry  Energy  Term:  The  correction  term  proportional  to 
( N-Z) 2/A  is  included  in  both  the  volume  energy  term  and  the 
surface  energy  term  in  order  that  "the  surface  and  volume 
energies  vary  in  a  parallel  manner  and,  in  particular, 
vanish  simultaneously  for  a  large  neutron  excess. 
(Conventinal  formulae  without  a  composition-dependent 
surface  energy  correspond  to  a  surface  tension  that  retains 
its  full  standard  value  even  when  the  neutron  excess  has 
made  the  volume  binding  vanish,  and  all  cohesion  has  been 
lost"). I*  "It  reflects  the  observation  that  for  a  given  A 
the  binding  energy  due  to  nuclear  forces  (i.e., 
disregarding  the  Coulomb  effect)  is  greatest  for  the 
nucleus  with  equal  number  of  neutrons  and  protons  and 
decreases  symmetrically  on  both  sides  of  N-Z.  The  simplest 
functional  form  expressing  these  empirical  facts  is  a  term 
in  (N-Z)2.  The  A-1  dependence  of  the  symmetry  energy  comes 
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about  because  the  binding  energy  contribution  per 
neutron-proton  pair  is  proportional  to  the  probability  of 
having  such  a  pair  within  a  certain  volume  (determined  by 
the  range  of  nuclear  forces),  and  this  probability  in  turn 
is  inversely  proportional  to  the  nuclear  volume.  The  extra 
stability  of  N=Z  nuclei  comes  about  at  least  in  part 
through  the  Pauli  exclusion  principle  since  two  identical 
nucleons  cannot  be  in  the  same  energy  state,  the  lowest 
state  for  a  given  number  of  nucleons  is  attained  (in  the 
absence  of  Coulomb  forces)  for  equal  numbers  of  neutrons 
and  protons."5 

(6)  Pai ring  Energy  Term:  This  correction  term  is  due  to  the 
coupling  of  pairs  of  identical  nucleons.  "This  term  Is  a 
quantitative  expression  of  the  fact  that  binding  energies 
for  a  given  A  depend  somewhat  on  whether  N  and  Z  are  even 
or  odd.  So  called  even-even  nuclei  (even  Z-even  N)  are  the 
stablest  and  for  them  the  pairing  energy  term  may  be  taken 
as  +11/A1/2;  for  even-odd  (even  Z-odd  N)  and  odd-even 
(odd  Z-even  N)  nuclei,  the  pairing  energy  term  Is  zero;  for 
odd-odd  nuclei  (odd  Z-odd  N)  the  pairing  energy  term  Is 
-11/A1/2.  The  greater  stability  of  nuclei  with  filled 
energy  states  is  apparent  not  only  in  the  larger  number  of 
even-even  nuclei  but  also  in  their  greater  abundance 
relative  to  the  other  types  of  nuclei.  On  the  average, 
elements  of  even  Z  are  much  more  abundant  than  those  of  odd 
Z  (by  a  factor  of  about  10).  For  elements  of  even  Z  the 
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isotopes  of  even  mass  (even  N}  account  in  general  for  about 
70  to  100  percent  of  the  element.  The  general  shape  of  the 
binding-energy  curve  (Fig.  2.3)  with  the  maximum  at  A  =  60 
comes  about  through  opposing  trends  with  mass  number  of  the 
relative  contributions  of  surface  energy  (decreasing  with 
A)  and  coulomb  and  symmetry  energies  (increasing  with 
A). "5 

The  semi -empi ri cal  mass  formula  can  be  written  as: 

Eb  >  Cj  A  [1  -  k(^)2]  -  c2  A2/3  [1  -  k (M)2] 

-  c3  Z2  A-1/3  +  c4  Z2  A-l  +  6  (2.11) 

"where  Eb  is  the  binding  energy,  that  is,  the  energy 
required  to  dissociate  the  nucleus  into  its  constituent 
nucleons,  and  A,  Z  and  N  are  the  mass  number,  proton  number 
and  neutron  number  respectively.  With  Eb  expressed  in  MeV 
the  coefficients  take  on  the  following  values: 
ci  =  15.677  MeV 
c2  *  18.56  MeV 
c3  *  0.717  MeV 
c4  =  1.211  MeV 
k  *  1.79 

6  *  pairing  term. 

The  equation  contains  only  six  empirically  adjusted 
parameters,  yet  this  simple  equation  yields  binding 
energies  that  agree  with  experimental  values  for  each  of 


the  approximately  1200  nuclides  of  known  mass  to  within 
less  than  10  MeV,  and  very  much  better  than  that  for  most. 
This  represents  a  remarkable  success  for  the  liquid-drop 
model  of  the  nucleus. "5 

The  mass  M  of  a  neutral  atom  whose  nucleus  contains  Z 
protons  and  N  neutrons  is: 

M  =  ZMU  +  N  M  -  E„  (2.12) 

H  n  d 

where  Eb  is  the  binding  energy,  Mh  is  the  mass  of  the 
hydrogen  atom  and  Mn  is  the  mass  of  the  neutron.  By 
determining  the  binding  energy  for  a  nuclide,  and  knowning 
its  constituents,  the  mass  may  easily  be  determined. 

2.4  Nuclear  Evaporation  Theory 

"The  application  of  quantum  mechanics  to  the  calculation  of 
nuclear  reactions  heavy  nuclei  gives  rise  to  extremely  complicated 
expressions  which  cannot  be  treated  by  ordinary  methods,  since  there 
is  no  approximate  solution  for  this  complex  many-bodied  problem  on 
account  of  the  Intense  interaction  between  the  constituents  of  atomic 
nuclei.  On  the  other  hand  Bohr  has  pointed  out  that  just  the  extreme 
facility  of  energy  exchange  gives  rise  to  a  characteristic 
simplification  of  the  course  of  every  nuclear  process  initiated  by  a 
collision  between  a  particle  or  a  light  quantum  and  a  nucleus.  It 
consists  in  the  possibility  of  dividing  the  process  into  two 
well -separated  stages.  The  first  is  the  formation  of  a  compound 
nucleus  in  a  well-defined  state  in  which  the  incident  energy  is 
shared  among  al1  the  constituents,  the  second  is  the  disintegration 


of  that  compound  system,  which  can  tc  treated  as  quite  independent  of 
the  first  stage  of  the  process.  This  conception  has  been  extremely 
fruitful  in  the  treatment  of  all  kinds  of  nuclear  reactions. 
Qualitative  statistical  conclusions  about  the  energy  exchange  between 
the  nuclear  constituents  in  the  compound  s^ii-e  have  led  to  simple 
explanations  of  many  characteristic  features  of  nuclear  reactions. 

In  particular,  the  use  of  thermodynamical  analogies  has  proved  very 
convenient  for  describing  the  general  trend  of  nuclear  processes. 

The  energy  stored  in  the  compound  nucleus  can  in  fact  be  compared 
with  the  heat  energy  of  a  solid  body  or  a  liquid,  and,  as  first 
emphasized  by  Frenkel,  the  subsequent  expulsion  of  particles  is 
analogous  to  an  evaporation  process. "12 

In  the  fission  process,  a  neutron  with  a  given  amount  of  kinetic 
energy  combines  with  a  heavy  nucleus  to  form  a  compound  system. 

After  the  formation  of  the  compound  nucleus  the  disintegration  stage 
may  be  treated  separately  in  most  cases  from  the  formation  stage  "in 
a  sense  that  the  mode  of  disintegration  of  the  compound  system 
depends  only  on  its  energy,  angular  momentum,  and  parity,  but  not  on 
the  specific  way  in  which  has  been  produced.  If  the  Incident 
particle  comes  within  the  range  of  forces,  its  energy  is  quickly 
shared  among  all  the  constituents  well  before  any  re-emission  can 
occur.  The  state  of  the  compound  system  Is  then  no  longer  dependent 
on  the  way  it  was  formed. "2 

Consider  the  reaction  a  +  X+  c+  Y  +  b,  where  c  is  the  compound 
nucleus.  For  the  evaporation  theory  to  be  used  to  express  the 
process  after  formation  of  the  compound  nucleus,  the  assumption  is 


made  that  "the  probability  that  the  compound  nucleus  c  decays  into 
channel  g  is  given  by: 


1  °c(3) 

G  (3)  =  — s -  where  (2.13) 

I  K2  0c(Y) 

Y 

6c(g)  =  Branching  ratio  (relative  decay  probability)  for 
the  compound  nucleus  c  to  decay  through  channel  g. 

Kg  =  iJigl  =  *g_1;  channel  wave  number  in  channel  g. 
oc(g)  =  Cross  section  for  the  formation  of  the  compound 
nucleus  through  channel  g. 

I  K  2oc(y)  =  Sum  of  K2oc  extended  over  all  channels  y. 

Y  Y 

into  which  c  can  decay. 

and  therefore  can  be  determined  from  the  cross  sections,  <jc."29 
The  energy  distribution  of  the  emitted  particles  is  first 
studied.  "The  kinetic  energy  in  channel  g  is  given  by: 

£B  *  c«  *  QaS  (2-14) 

a  =  entrance  channel 


g  =  exit  channel 


Qag  =  Q  value  for  the  reaction  *  Ea*  +  Ea"  -  Eg*  -  Eg"  (single  ' 
indicates  states  of  nuclei  and  double  '  indicates  the  states 


of  the  Incident  or  emerging  particle.) 

Since  most  of  this  energy  is  the  kinetic  energy  of  the  emitted 
particle  b,  and  only  a  small  part  is  the  recoil  energy  of  the 
residual  nucleus  Y,  we  consider  eg  the  kinetic  energy  of  b. 
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We  write  here:  eg  =  egy  -  Eg*,  where  egy  =  e  +  Qab  and  e  =  ea0 
is  the  kinetic  energy  in  the  entrance  channel  (both  a  and  X  in  their 
ground  states),  and  represents  the  maximum  value  of  eg,  which  occurs 
if  the  residual  nucleus  Y  remains  in  its  ground  state.  Eg*  is  the 
excitation  energy  of  Y  after  the  reaction  and  the  difference  between 
the  energy  of  Y  after  the  reaction  and  the  energy  of  the  ground  state 
of  Y.  To  every  level  Eg*  of  the  residual  nucleus  there  corresponds 
an  energy  eg  of  the  outgoing  particle,  as  long  as  Eg*  <  egy."13 

The  formation  cross  section  oc  generally  increases  with  in¬ 
creasing  channel  energy  (owing  to  better  penetration  of  the  barrier), 
except  for  neutrons  with  i  =  0  where  oc  decreases,  but  very  slowly 
(proportionally  to  k_1  or  e-1/2).  In  all  cases  the  product  k2oc 
which  appears  in  Eq.  (2.13)  is  a  monotonical ly  increasing  function  of 
the  channel  energy  e.  If  the  energy  is  high  enough,  a  great  many 
levels  of  the  residual  nucleus  can  be  excited.  The  enerq*  distribu¬ 
tion  of  the  particles  becomes  continuous  when  the  energies  Eg  !  are 
closer  together  than  the  definition  of  the  energy  in  the  incident 
beam,  or  closer  than  the  energy  resolution  of  the  experimental  ar¬ 
rangement  for  detection  of  the  reaction  products.  The  shape  of  the 
distribution  function  Gb(e)de  for  the  number  of  particles  b  emitted 
between  e  and  e  +  de  is  then  given  by: 

6bU)de  *  l  Gc(a)  (2.15) 

e  <  e3  <  e  +  de 

where  the  sum  is  extended  over  channels  g  whose  energy  lies  within 
the  energy  interval  de.  The  number  of  terms  in  this  sum  is  given  by 


V.  •/,  *_• 


the  number  of  levels  of  the  residual  nucleus  Y  with  an  excitation 
energy  Eg*  between  E  and  E  -  de,  where  E  =  eby  -  e.  We  call  this 
number  Wy(E)de,  and  Wy(E)  the  level  density.^ 

The  level  density  can  be  used  to  determine  the  relative 
intensity  distribution  of  the  outgoing  particles  Ib(e)de  by  the 
rel ationship: 

Ib(e)de  =  constant  e  a  (8)  WyUbY  -  e)  dc  (2.16) 

wnere  ac(6)  =  ocU)  is  a  function  of  the  channel  energy  c  =  eB. 

"We  obtain  a  rough  estimate  of  the  shape  of  Ib(e),  especially 
for  neutrons,  by  introducing  the  following  Taylor  expansion  of  the 
logarithm  of  the  level  density: 

n(E)  =  log  W(E)  (2.17) 

around  the  maximum  energy  eby  by  which  the  residual  nucleus  Y  can  be 
excited. 

v  UbY  -  e)  *  V  (eb.  )  -  e  (-gjj-)  (E«ebY)  +  .  .  .  (2.18) 

If  the  expansion  Is  used  to  approximate  W(egy  -  e)  In  the  relation¬ 
ship  for  Ib(e)de,  we  can  abosrb  the  factor  coming  from  n  UbY  -  e) 
Into  the  constant  and  get: 


Ib  U)de 


constant  e  ar(e)  exp  -  -r-A — r  dc 

C  9  \EbY; 


where  the  function  q  is  determined  by: 

1  dn 

0*dE 


(2.19) 


(2.20) 


0  has  the  dimensions  of  an  energy  and  can  be  interpreted  as  a  nuclear 
temperature.  This  is  the  well-known  thermodynamical  relation  between 
entropy  U)  and  temperature  (e)."^ 


The  application  of  the  temperature  concept  to  nuclear  reactions 


can  be  understood  in  the  following  way:  The  incident  particle  'a' 
enters  the  target  nucleus  and  forms  a  highly  excited  compound 
nucleus.  The  excitation  energy  can  be  considered  heat  energy 
delivered  by  the  impact  of  'a'  on  the  target.  The  heating  of  the 
compound  system  causes  an  evaporation  of  neutrons  or  other  particles 
and  the  energy  distribution  of  theemitted  neutrons  is  similar  to  the 
Maxwell  Distribution.  The  potential  barrier  distorts  the  Maxwell 
distribution  for  charged  particles  by  depressing  the  emission  of  low 
energy  particles. 

The  temperature  which  determines  the  Maxwell  distribution  is 
given  by  0  Uby): 

0  (ebYJ  *  (dn/dt) :  E  *  cbY  (2,21 

It  is  the  temperature  of  the  residual  nucleus  Y  at  the  excita¬ 
tion  energy  ebY  (i.e.,  the  temperature  of  the  residual  nucleus  after 
emission).  In  the  case  of  the  fissioning  of  U-236  into  two  frag¬ 
ments,  the  excitation  energy  of  the  fragments  will  be  the  important 
factor  in  determining  the  energy  allowed  for  neutron  emission. 


CHAPTER  3 


CODE  DEVELOPMENT 

The  code  to  determine  the  prompt  neutron  spectra  and  neutron 
multiplicity  for  the  thermal  fission  of  U-235  is  comprised  of  a  main 
program,  a  block  data  subprogram,  a  mass  correction  subroutine,  a 
binding  energy  subroutine,  and  a  neutron  kinetic  energy  subroutine. 

3.1  Block  Data  Subprogram 

The  Block  Data  Subprogram  is  used  ti  initialize  the  large 
arrays,  AN(IA)  ZP(IA),  YA(IA),  and  ZA(IA),  which  are  the  95  values  of 
the  atomic  number  of  the  IA-th  nuclide  chain,  the  most  probable 
charge  number  of  the  IA-th  nuclide  chain,  the  chain  yield  of  IA-th 
nuclide  chain,  and  the  starting  proton  number  for  the  IA-th  chain. 

The  values  for  AN(IA)  range  from  71.0  to  165.0,  which  is  the  atomic 
number  range  for  the  fission  fragments  of  U-235.  The  values  of 
ZP(IA)  were  taken  from  Table  52.10,  Reference  [7],  the  values  for 
YA(IA)  were  found  in  Reference  [14],  The  values  of  ZA(IA)  were  taken 
from  the  chart  of  the  nuclides  by  examining  each  nuclide  chain  and 
extracting  the  stable  species  of  highest  proton  number  (lowest 
neutron  number).  These  values  are  passed  to  the  main  program  in  a 
common  statement  called  COMMON/BLK/. 

3.2  The  Average  Heavy  and  Light  Masses 

The  charge  number  of  the  IZ-th  nuclide  with  mass  number  AN(IA) 
can  be  placed  into  an  array  ZN(IA,IZ)  by  using  two  DO-LOOPS  over  IA 
(number  of  nuclide  chains)  and  IZ  (number  of  nuclides  in  each  chain) 


The  weighted  averages  are  then  printed  as  output,  to  give  the  average 
heavy  fragment  mass  and  the  average  light  fragment  mass. 


3.3  Determination  of  Fragment  Kinetic  Energy  Constant 

The  formula  for  the  electrostatic  repulsion  of  two  charged 
spheres  is  used,  with  minor  modification,  to  approximate  its  fragment 
kinetic  energy  (EF). 


EFtotal 


Z(l)  *  Z ( 2 )  *  e2 
RR[(A(l)i/3  +  (A{2))i/3] 


(3.3) 


where  Z(l),  A(l)  are  the  proton  number  and  atomic  mass  number  of 
fragment  1  and  Z(2),  A(2)  are  the  proton  number  and  atomic  mass  num¬ 
ber  for  fragment  2,  e  is  the  elementary  charge,  and  RR  is  the  effec¬ 
tive  distance  (see  Fig.  3.1)  between  the  centers  of  the  fragments  at 
the  moment  the  neck  of  the  fissioning  system  splits  and  coulomb 
repulsion  forces  the  fragments  to  repel  away  from  each  other.  The 
value  for  RR  is  determined  by  plotting  results  for  different  values 
and  obtaining  an  empirical  "best  fit".  The  value  for  e  is 
4.80298  x  1010  electrostatic  units.  The  conversion  constant  is  used 
to  transform  ergs  to  MeV  1.60219  x  10"6  ergs/MeV.  Therefore  we  can 
state: 


EFtotal 


Z(l)  *  Z ( 2 ) 

RRC(A(l))i/3  ♦  (A(2) ) l/3] 


*  Constant 


(3.4) 


where  the  constant  is  equal  to  1.439899  and  RR  is  in  Fermis  (10" i 3 
centimeters) . 

3.4  The  Main  Program  Loop 

The  greatest  portion  of  the  main  program  consists  of  two 
D0-L00PS,  one  over  IA  (number  of  nuclide  chains)  and  IZ  (number  of 


(a)  Uranium  235  +  neutron  -►  Uranium  236 


(b)  Upon  separation  during  fission  process  a  "neck  Is  formed" 


(c)  The  two  fragments  are  separated  by  an  "effective  distance 
RR  determined  In  this  program  to  be  1.835  fermls  (1.835  x 
10" 13  centimeters) 


Fig.  3.1  Effective  Distance  Between  Fission  Fragments 


nuclides  in  the  chain),  which  in  effect  performs  all  operations  in 
the  program  for  each  possible  fission  fragment  released  in  the  binary 
fission  (thermal  neutron  induced)  of  Uranium-235.  The  index  IA  is 
transformed  to  ID  =  NA  -  IA  +  1  (where  NA  =  total  number  of  chains) 
so  that  the  heaviest  chain  Is  used  first,  and  IZ  is  transformed  to 
IE  =  NZ  -  IZ  +  1  (where  NZ  =  total  number  of  nuclides  in  the  chain), 
so  the  smallest  proton  number  is  used  first,  and  hence  the  largest 
neutron  number  nuclide  is  used  to  facilitate  the  neutron  emission 
scheme  used  in  the  code. 

The  atomic  number  of  fragment  1  is  stored  in  A(l)  and  is  equal 
to  AN ( ID) .  The  atomic  number  of  fragment  2  is  stored  in  A(2)  and  is 
equal  to  236.0  -  AN(ID)  so  that  the  two  fragments  add  up  to  the  com¬ 
plex  nucleus  U-236's  total  number  of  nucleons  which  is  236.  The  pro¬ 
ton  number  of  fragment  1  is  stored  in  Z(l)  and  is  equal  to  ZN( ID,IE) 
and  the  proton  number  of  fragment  2  Is  equal  to  92.0  -  ZN(ID,IE)  so 
that  the  sum  of  the  proton  numbers  Is  equal  to  92.0  which  is  the  pro¬ 
ton  number  for  the  complex  nucleus  U-236.  The  binding  energy  sub¬ 
routine  is  called  and  the  rest  mass  for  fragment  1  and  fragment  2  is 
returned  so  that  the  energy  released  ( ER )  in  fission  can  be  deter¬ 
mined  by  subtracting  the  rest  mass  of  the  two  fragments  from  the  rest 
mass  of  the  complex  nucleus.  By  using  Z(l),  Z(2),  A(l),  A(2),  the 
fragment  kinetic  energy  can  be  determined  for  the  two  fragments  under 
consideration.  In  the  main  program,  all  ER  and  all  EF  are  summed  and 
an  average  ER  and  average  EF  can  be  computed  after  completion  of  the 
DO  LOOP.  The  incident  neutron,  being  of  thermal  energy,  adds  .025  eV 
to  the  system.  Although  this  is  an  insignificant  amount,  it  is  in- 


eluded.  The  total  excitation  energy  (EE( ID, IE) )  of  the  fragments  can 
be  determined  by: 

EE(ID.IE)  =  ER  +  EN  -  EF  (3.5) 

where  EN  is  the  incident  neutron  energy.  The  values  for  EE(ID,IE) 
are  also  summed  up  in  order  to  determine  an  average  excitation  energy. 

The  binding  energy  subroutine  is  called  a  second  time  and  enters 
the  subroutine  with  the  atomic  mass  numbers  and  proton  numbers  of 
fragment  1  and  fragment  2.  For  each  of  these  fragments,  the  rest  mass 
is  determined,  and  the  rest  mass  RM(I),  neutron  separation  energy  SEN, 
SEN2  and  pairing  energy  T5(I)  is  determined  for  9  successive  neutron 
emissions  in  order  to  form  neutron  decay  chains  which  can  be  examined 
to  determine  the  number  of  neutrons  emitted  within  the  allowable 
energy  of  the  system. 

Fragment  1  (the  heavy  fragment)  is  examined  first.  The  excita¬ 
tion  energy  on  the  heavy  fragment  (EEH)  can  be  computed  by  multiplying 
the  excitation  energy  of  the  system  by  the  fragment  atomic  mass  number 
and  dividing  by  the  atomic  mass  number  of  U-236. 

EEH  *  EE( ID, IE )*A( 1 )/236 .00  (3.6) 

Here,  it  is  assumed  tht  all  236  nucleons  are  in  thermodynamic  equi¬ 
librium  and  on  the  average,  each  nucleon  shares  an  equal  amount  of 
residual  energy  in  the  nucleus.  The  separation  energy  of  a  neutron 
(SEN)  can  be  computed  by: 

SEN  -  939.553  -  (RM( I)  -  RM(I  +  1))  (3.7) 

where  939.553  is  the  rest  mass  of  the  neutron,  RM(I)  is  the  rest  mass 
of  the  fragment  and  RM(I  +  1)  is  the  rest  mass  of  the  next  fragment  in 
the  neutron  decay  chain. 


The  maximum  kinetic  energy  (EKMAX)  of  the  emitted  neutron  from 


the  heavy  fragment  is  the  excitation  energy  of  the  heavy  fragment 
minus  the  separation  energy  of  the  neutron  minus  the  pairing  energy: 

EKMAX  =  EEH  -  SEN  -  T5(I)  (3.8) 

The  energy  subroutine  is  called  to  determine  a  neutron  energy, 
based  on  evaporation  theory,  for  the  emitted  neutron.  If  the  neutron 
is  emitted,  the  multiplicity  NU( ID, IE)  for  the  fragment  is  increased 
by  one.  The  excitation  energy  for  the  fragment  (EEH)  is  now  equal  to 
the  maximum  kinetic  energy  of  the  emitted  neutron  minus  the  kinetic 
energy  (EK( ID, IE, I ) )  determined  from  the  energy  subroutine.  This 
process  is  continued  until  the  excitation  energy  (EEH)  is  less  than 
the  sum  of  the  separation  energy  of  the  neutron  and  the  pairing  energy 
of  the  fragment.  When  this  occurs,  the  residual  excitation  energy 
(EEH)  can  be  set  equal  to  the  energy  available  for  Gamma  Decay  (GEH) 
from  the  heavy  fragment. 

The  complementary  (light)  fragment  also  Is  considered  in  the  same 
manner  as  the  heavy  fragment.  In  this  portion  of  the  program,  the 
excitation  energy  of  the  light  fragment  is  EEL,  and  the  separation 
energy  is  SEN2. 

The  index  counter  (I)  ranges  from  1  to  20.  The  numbers  1  through 

10  are  used  for  the  heavy  fragment  decay  chain  and  number  20  down  to 

11  are  used  for  the  light  fragment  decay  chain. 

The  Gamma  decay  energy  for  the  light  fragment  Is  called  GEL.  The 
total  energy  available  for  gamma  emission  (GE)  is  then: 

GE  *  GEH  +  GEL  (3.9) 


Vv, 


The  neutron  multiplicity  can  be  multiplied  by  the  fractional 
yield  of  the  nuclides  and  summed  over  all  nuclides  to  get  a  yield 


weighted  neutron  multiplicity.  Also,  a  yield  weighted  gamma  energy  is 
obtained. 

This  completes  the  main  loop,  and  the  values  of  the  average 
release  energy  (ER),  fragment  total  kinetic  energy  (EF),  excitation 
energy  (E),  Gamma  energy  (GE)  and  the  yield  weighted  neutron 
multiplicity  is  sent  to  output. 

In  order  to  produce  a  neutron  kinetic  energy  spectrum,  20  energy 
bins  are  established  in  0.5  MeV  increments.  The  number  of  neutrons  in 
each  of  these  bins  are  stored  in  NN( 1 )  through  NN( 20) .  This  data 
forms  the  basis  for  a  spectrum  of  the  number  of  neutrons  per  energy 
bin  versus  the  energy  of  the  bins. 

3.5  Binding  Energy  Subroutine 

The  binding  energy  subroutine  utilizes  the  Meyer-Swlatecki  mass 
formula  to  compute  the  rest  mass  and  binding  energies  for  the  fission 
fragments  in  order  to  determine  the  release  energy  of  the  fission 
process.  The  rest  mass  of  the  fragments  and  the  successive  daughters 
of  the  neutron  evaporation  decay  chain  Is  computed  in  order  to 
determine  the  separation  energy  of  the  emitted  neutron. 

The  subroutine  is  entered  with  the  atomic  number  A(I)  of  the 
heavy  (1*1)  fragment  and  light  (1*2)  fragment,  Z(I)  of  the  heavy 
fragment  and  light  fragment,  from  which  the  neutron  number  U(I)  is 
determined  (U(I)  *  A(I)  Z(I)). 


The  first  time  this  routine  is  called,  counter  J  is  set  equal  to 
2,  denoting  that  the  loop  for  rest  mass  calculation  will  be  run  twice, 
once  for  the  heavy  fragment,  once  for  the  light  fragment. 

The  rest  mass  loop  consists  of  computing  the  volume  energy  term 
Tl,  the  surface  energy  term  (T2) ,  the  coulomb  energy  term  (T3),  the 
coulomb  correction  term  (T4)  and  the  pairing  energy  term  T5(I).  The 
pairing  energy  term  is  placed  in  a  vector  to  be  used  in  the  main 
program  for  determining  the  maximum  kinetic  energy  of  the  emitted 
neutron  (EKMAX  =  EEH  -  SEN  -  T5( I ) ) .  The  pairing  term  can  be  one  of 
three  different  relationships  depending  on  whether  A(I)  is  odd  or  even 
and  whether  Z(I)  is  odd  or  even,  so  a  test  is  performed  to  determine 
whether  A  and  Z  are  odd  or  even  and  the  appropriate  pairing  term  rela¬ 
tionship  is  used.  The  Binding  energy  (BE)  is  determined  by  summing 
the  five  terms  and  the  units  are  in  MeV.  In  order  to  determine  the 
binding  energy  in  amu  (BEA)  the  binding  energy  Is  divided  by  931.5016 
MeV  per  amu.  The  rest  mass  in  MeV  Is  equal  to  the  negative  of  the 
binding  energy  in  MeV  plus  the  proton  number  times  the  rest  mass  of 
the  proton  plus  the  neutron  number  times  the  rest  mass  of  the  neutron. 

RM( I)  *  BE( I)  +  Z( 1)  *  938.7906  +  (A(I)  -  Z( I) )  *  939.553  (3.10) 

The  rest  mass  in  amu  Is  determined  by  dividing  the  rest  mass  in  MeV  by 
931.5016  MeV  per  amu.  In  order  to  apply  the  rest  mass  corrections  in 
MeV  for  proton  and  neutron  number  as  determined  in  the  rest  mass  cor¬ 
rection  subroutine  (MASCOR)  the  terms  YB(IMC)  (neutron  mass  correc¬ 
tion)  and  YD(IMA)  (proton  mass  correction)  are  added  to  the  rest  mass 
in  MeV  ( RM( I ) ) .  The  subroutine  returns  to  the  main  program  (after  the 


first  call)  with  the  rest  masses  and  pairing  energies  of  the  heavy  and 
light  fragments. 


The  second  call  of  the  binding  energy  subroutine  enters  with  the 
atomic  number  AN(ID)  and  proton  number  ZN(ID,IE)  of  the  heavy  frag¬ 
ment.  The  neutron  number  U(I)  is  determined  by  subtracting  the  proton 
number  from  the  atomic  number.  The  heavy  fragment  and  its  nine 
successive  evaporation  decay  daughters  are  stored  in  A(l)  through 
A ( 10 )  with  corresponding  Z(l)  through  Z( 10)  (proton  numbers)  and  U ( 1 ) 
through  U( 10 )  (neutron  numbers).  The  complementary  light  fragment  is 
stored  in  A(20)  with  a  corresponding  Z( 20)  and  U(20).  The  nine 
successive  decay  products  are  stored  in  A( 19)  through  A(ll)  with 
corresponding  Z( 19 )  through  Z(ll)  (proton  numbers)  and  U( 19)  through 
U ( 1 1 )  (neutron  numbers).  The  rest  mass  is  then  determined  for  each  of 
these  fragments,  requiring  a  20  count  loop,  and  this  Information  Is 
passed  back  into  the  main  program  so  that  the  separation  energy  for 
the  heavy  fragment  (SEN  =•  939.553  -  (RM(I)  -  RM(I  +  1))  and  the 
separation  energy  for  the  light  fragment  (SEN2  *  939.553  -  (RM(I)  - 
RM(I  -  1))  can  be  determined  for  use  In  the  primary  loop  of  the  main 
program. 

3.6  Energy  Subroutine 

The  energy  subroutine  uses  an  evaporation  theory  model  of  deter¬ 
mining  the  kinetic  energy  of  the  emitted  nucleon.  This  subroutine  Is 
based  on  Its  work  of  I.  Dostrovsky,  Z.  Fraenkel  and  6.  Friedlander, 
which  Is  based  on  the  model  of  Weisskopf,  and  uses  a  Monte  Carlo 
selection  of  the  kinetic  energy  based  on  the  evaporation  spectrum. 


The  required  input  to  the  subroutine  is  the  Atomic  Number  of  the 


fission  fragment  Al(II),  the  proton  number  of  the  fission  fragment, 
Zl(II),  the  neutron  number  of  the  fission  fragment  Ul(II),  the 
excitation  energy  of  the  fission  fragment  (EEH  for  the  heavy  fragment 
and  EEL  for  the  light  fragment),  and  the  separation  energy  of  the 
fission  fragment  (SEN  for  the  heavy  fragment,  and  SEN2  for  the  light 
fragment).  The  subroutine  initially  determines  the  level  density 
parameter  AU  by  the  following  empirical  relationship: 

AU  =  A( (1.  -  1.3  *  THETA/A1 ( I ) )**2)  (3.11) 

where  A  =  Al(II)/20.0,  the  atomic  number  divided  by  20.0.  A  is  the 
uncorrected  level  density  parameter.  THETA  is  the  neutron  number 
minus  the  proton  number  divided  by  the  atomic  number  (THETA  *  (Ul(II) 

-  Zl( II) )/Al ( II) ) .  A  test  is  performed  to  determine  if  Al(II)  and/or 
Zl(II)  is  odd  or  even  so  that  the  appropriate  pairing  term  DELTA  can 
be  determined. 

If  the  subscript  II  <  10,  this  denotes  the  heavy  fragment,  and  if 
II  >  11,  this  denotes  the  complementary  light  fragment.  Initially, 
the  heavy  fragment  Is  used  and  the  maximum  kinetic  energy  RJ1  for  the 
neutron  emitted  Is  computed  by: 

RJ1  =  EEH  -  SEN  -  DELTA  (3.12) 

The  neutron  kinetic  energy  selection  scheme  uses  the  substitution 
XX  *  EK(ID2,  IE2,  II)  +  BETA  where  EK  is  the  kinetic  energy  of  the 
emitted  neutron  and  BETA  is  an  empirical  formula 

BETA  *  (2.12  *  A  **  (-2/3)  -  0.050)/(0.76  +  2.2  *  A  **  (-1/3)) 


With  the  inverse  cross  section  as  given  by  oc/ag  =  a(l  +  3e) 
where  a  =  0.76  +  2.2  *  A  **  (-1/3)  and  e  given  as  BETA  above,  and  <jc 
and  Og  are  the  capture  and  geometric  sections  respectively,  and  the 
level  density  is  given  by: 

W(E)  =  C*EXP  (2.  *(AU(EE  -  DELTA))  8**  1/2) 
where  DELTA  =  0  for  odd-odd  nuclei,  and  DELTA  *  0  for  other  types. 

The  following  equation  is  obtained  for  neutron  emission, 
g  m  r2A  2/3 

PX( e)de  =  - -  exp  {-2*[AU  (EE  -  DELTA  ))**  i}* 

TT  h2/3  0  0  2 

e*a  (1  +  ^~)  *  exp  {2  *  [AUn  *  (EE  -  SEN  -  DELTAn  -  e ) ]**  }}de 

(3.14) 

where  the  subscript  o  and  n  refer  to  the  original  and  residual  nucleus 
respectively  and  e  is  the  kinetic  energy  of  the  neutron  emitted.  To 
simplify  this  notation,  the  substitution  XX  *  e  +  BETA  Is  made.  By 
differentiation,  the  value  XMAX  for  which  PX(XX)  Is  a  maximum  may  be 
found  to  be: 

XMAX  *  [(AU  *  RJ1  +  0.25)  **  0.5  -  0.5J/AU  (3.15) 

and  now  PX(XX)  may  be  normalized  to  PX(XMAX)  *  1;  it  then  becomes: 

PX(XX)  »  «p  {AU  *  XMAX  -  l(AU  *  (RJ1  -  XX))  **  1/2J}"16 

(3.16) 

"The  actual  choice  of  a  particular  value  of  XX  Involves  the 
consecutive  drawing  of  two  numbers  between  0  and  1.  The  first  random 
number  RND(l)  is  used  to  choose  a  value  of  XX  between  0  and  RJ1  (the 
maximum  possible  kinetic  energy  for  the  given  particle: 


XX  =  RN0(1)  *  RJ1 


(3.17) 
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The  second  random  number  RND(2)  determines  whether  the  value  of  XX  is 
to  be  accepted.  If  PX(XX)  >  RND(2),  accept  the  value  of  XX.  If 
PX(XX)  <  RND(2),  reject  the  value  of  XX  and  select  another  set  of 
random  numbers. "15 

If  XX  -  BETA  is  negative,  another  set  of  random  numbers  is 
selected  to  prevent  negative  energies  which  are  not  physically 
meani ngful . 

"The  range  of  possible  kinetic  energies  of  the  outgoing  particle 
was  divided  into  three  regions,  and  a  different  procedure  was  used  for 
each,  (a)  For  the  region  XX  <  2  *  XMAX,  the  normalized  PX(XX)  was 
computed  from  Eq.  (3.16)  and  the  value  compared  with  RND(2).  If  the 
computed  value  was  smaller  than  RND(2),  it  was  rejected  and  a  new 
random  number  set  (RND(l),  RND(2))  is  chosen.  If  the  computed  value 
of  PX(XX)  was  greater  than  RND(2) ,  the  kinetic  energy  chosen  by  the 
first  random  number  RND(l)  was  accepted  and  its  value  computed: 

EK  *  XX  -  BETA  (3.18) 

(b)  In  the  second  region,  2  *  XMAX  <  XX  <  15  *  XX,  use  was  made  of 
the  fact  that,  for  any  positive  value  of  XX,  the  probability  given  by 
a  Maxwell  distribution: 

XX 

PX(XX)  *  x^x  exp  [1  -  (XX/XMAX)]  (3.19) 

is  greater  than  (or  for  XX  *  XMAX  equal  to)  that  given  by 

Eq.  (3. 16). "15  for  this  region  the  same  procedure  used  in  (a)  above 

was  used  but  PX(XX)  was  calculated  using  Eq.  (3.19). 


(c)  For  the  region  XX  >  15  *  XMAX,  values  were  all  rejected 
because  of  the  very  small  probability  of  neutrons  being  found  in  this 
region. 

Once  the  value  of  XX  is  accepted,  then  the  kinetic  energy  EK(ID2, 
IE2,  II)  is  calculated  by  EK ( I D2 ,  IE2,  II)  =  XX  -  BETA. 

The  value  for  EK ( ID2 ,  IE2,  II)  is  in  the  center  of  mass 
coordinate  system  and  must  be  transformed  into  the  laboratory  system. 

In  the  main  program,  the  correction  factors  (for  the  heavy 
fragment,  CL1  and  for  the  light  fragment  CL2)  are  computed  and  passed 
into  the  subroutine  by  a  common  statement. 


m  _  939.553  *  RM(2)  .  rr 

CL1  RM(l)  (RM(l)  +  RM(2))  CF 

(3.20) 

r|„  939.553  *  RM(1) 

CL2  RM(2)  (RM(1)  +  RM(2) )  EF 

(3.21) 

where  RM(I)  is  the  rest  mass  of  the  heavy  (1)  or  the  light  (2) 
fragment  and  EF  is  the  total  kinetic  energy  of  the  fragments  and 
939.553  Is  the  rest  mass  of  the  neutron. 

The  laboratory  kinetic  energy  of  the  neutron  being  emitted  is 

then: 

EK ( ID2,  IE2,  II)  -  CL1  ♦  EK( ID,  IE,  II)  ♦  2.0  *  SQRT(CL1  *  EK( ID2, 
IE2,  II))  *  COS(ANGLE)  (3.22) 

where  ANGLE  is  a  randomly  distributed  angle  between  0  and  2w: 

ANGLE  =  RND(3)  *  6.283185  .  (3.23) 

The  light  fragment  also  has  a  laboratory  system  kinetic  energy: 
EK( ID2,  IE2,  II)  *  CL2  ♦  EK( ID,  IE,  II)  ♦  2.0  *  SQRT(CL2  * 

EK( ID2,  IE2,  II))  *  C0S( ANGLE) 


.w.vyv/.vv 


Once  these  neutron  kinetic  energies  are  determined,  a  common 
statement  passes  them  back  into  the  main  program  and  the  new  residual 
energy  (EEH,  EEL)  is  found. 

3.7  Subroutine  for  Mass  Corrections 

The  MASCOR  subroutine  performs  a  cubic  spline  fit  of  selected 
data  points  from  the  mass  correction  graphs  for  proton  and  neutron 
number  as  plotted  by  Myers  and  Swiatecki.  The  graph  used  for  this 
subroutine  is  the  experimental  shell  effect  graph,  Mass 
(experimental )-Mass  (liquid  drop),  and  the  discrepancy  between  the  two 
values  is  in  MeV).  The  routine  is  performed  for  both  neutron  number 
and  proton  number  and  these  correction  terms  are  then  transferred  to 
to  main  program  and  back  into  the  binding  energy  subroutine  where  the 
values  are  added  to  the  rest  mass  (in  MeV)  terms. 

If  the  set  of  known  points  in  Increasing  order  are  (UX(1), 

YM(1) ) ,  (UX(2) ,  YM ( 2 ) ) . (UX(NN),  YM(NN)),  then  a  spline  fit  is 

accomplished  by  connecting  each  pair  of  adjacent  points  with  a  section 
of  a  third-degree  polynomial,  matching  up  the  sections  so  that  the 
first  and  second  derivatives  are  continuous  at  each  point.  Let  ZK(1), 

ZK(2) . ZK(NN)  be  the  values  of  the  second  derivative  at  the 

points.  Then  between  points  (UX(K),  YM(K))  and  (UX,K+1),  YM(K+1)), 
the  second  derivative  has  the  value 

YM"  «  ZK(K)  U*  +  ZK(K+1)  UX 

where  DK  »  UX(K*1)  -  UX(K)"16 


(3.25) 


Integrating  Eq.  (3.24)  will  result  in  getting  a  first  derivative: 


YM'  =  -  ZK(K)  [(UX(K+1)  -  UX)2/2DK]  +  ZK(K+1)  [(UX  -  UX(K))2/2DK] 
+  CN(1)  (3.26) 

where  CN(1)  is  the  constant  of  integration. 

Integrating  once  more  will  result  in  the  equation  of  the  curve: 

YM  =  ZK(K)  [(UX(K-l)  -  UX)3)/6DK]  +  ZK(K+1)  [(UX  -  UX(K))3/6DK] 

+  CN(1)  *  UX  +  CN( 2)  (3.27) 

and  again  CN(2)  is  the  constant  obtained  by  integration. 

Knowing  that  the  curve  passes  through  (UX(K),  YM(K))  and 
(UX(K+1),  YM(K+1)),  the  values  of  the  integration  constants  may  be 
obtained. 

YM(K)  =  (ZK(K)  *  (DK  **  2)/6)  +  CN(1)  *  UaU)  +  CN(2)  (3.28) 

YM(K+1 )  *  (ZK(K+1)  *  (DK**20/6)  +  CN(1)  *  UX(1)  +  CN(2)  (3.29) 

from  which  the  values  of  CN( 1 )  and  CN(2)  are: 

CN( 1 )  =  [(YM(K+1)  -  YM(K) )/DK]  -  [ZK(K+1)  -  ZK(K))  *  DK/6] 

(3.30) 

CN(2)  *  [( YM(K)  *  UX(K+1)  -  YM(K+1)  *  UX))/DK]  - 

[ ( ZK( K)  *  UX(K+1)  -  ZK(K+1)  *  UX(K))  *  DK/6]  (3.31) 

If  the  values  for  CN(1)  and  CN(2)  are  substituted  into  Eq.  (3.27) 
the  following  expression  is  obtained 

YM  *  [ZK(K)(UX(K+1)  -  UX)3/6DK]  +  [ZK(K+1)(UX  -  UX(K))3/6DK] 

♦  L(UX(K+1)  -  UX)(YM(K)/DK  -  ZK(K)  *  DK/6) J 
+  L(UX  -  UX(K) )  ( YM(K+1)/DK  -  ZK(K+1)  *  DK/6)]  (3.32) 


All  the  values  in  Eq.  (3.32)  are  known  with  the  exception  of  ZK(K)  and 
ZK(K+1),  which  are  the  second  derivatives  at  the  endpoints  of  the 
interval . 


One  condition  which  will  help  determine  these  values  is  that  the 
slope  at  (UX(K),  YM(K))  as  determined  from  Eq.  (3.26)  must  be  the 
same  as  that  determined  by  the  corresponding  formula  for  the  interval 
(UX(K-l) ,  YM(K-l) )  to  (UX(K) ,  YM(K)).  When  the  value  of  CN(1)  from 
Eq.  (3.30)  is  used  in  Eq.  (3.26)  the  equation  becomes:"16 

YM*  =  -  [ ZK ( K ) ( UX ( K+ 1 )  -  UX)2/2DK]  +  [ZK(K+1)(UX  -  UX(K) )2/2DK] 

+  [( YM(K+1)  -  YM(K) )/DK]  -  [(ZK(K+1)  -  ZK(K))DK/6]  (3.33) 

For  the  interval  preceding  this  one: 

YM'  =  -  [ZK( K— 1 ) ( UX( K)  -  UX ) 2/2DK( K-l ) ]  +  [ZK( UX 

-  UX(K-l) )2/2DK(K-l)]  +  [( YM(K)  -  YM(K-1)/DK(K-1)] 

-  [( ZK(K)  -  ZK(K-l) )  *  DK(K-l)/6]  (3.34) 

At  the  point  (UX(K),  YM ( K ) )  the  following  relationship  Is  derived  by 
setting  Eq.  (3.33)  equal  to  Eq.  (3.34). 

ZK(K-l) (DK(K-l)/6)  +  ZK[(DK(K-1)  +  DK)/3]  +  ZK(K+1)  •  (DK(K)/6) 

»  C(YM(K+1)  -YM(K) )/DK]  -[( YM(K)  -YM(K-l) )/DK(K-l)] 

(3.35) 

For  each  of  the  interval  points,  K  *  2,3,...,NN-1  there  is  an 
equation  like  Eq.  (3.35).  "There  are  NN-2  equations  in  the  NN 
unknowns  ZK(1),  ZK(2),  ...  ZK(NN).  Two  or  more  conditions  may  be 
specified  in  order  to  determine  these  equations  completely.  It  is 
customary  to  place  some  additional  condition  on  ZK ( 1 )  and  ZK(NN),  the 
values  of  the  second  derivative  at  the  endpoints.  There  are  several 
reasonable  choices  for  these  values,  and  the  particular  choice  will 
influence  the  shape  of  the  fit,  especially  near  the  endpoints.  We 
will  require  the  third  derivative  to  be  continuous  at  (UX(2),  YM(2)) 


YM*  "  =  -  ZK(K)/DK{K)  +  ZK(K+1)/DK(K)  (3.36) 

If  the  values  for  K  =  1  and  K  =  2  are  equated,  as  well  as  for  K  =  NN-2 
and  K  =  NN-1,  we  obtain: 

-  ZK ( 1 ) / DL ( 1 )  +  ZK ( 2 ) ( 1 / OK ( 1 )  +  1/DK(2) )  -  ZK(3)/DK(2)  =  0 

(3.37) 

and 

-  ZK( NN-2) /DK( NN-2)  +  ZK(NN-l) (l/DK(NN-2)  +  1/DK(NN-1)) 

-  ZK ( NN ) / DK ( NN- 1 )  =  0  (3.38) 

The  above  two  equations  along  with  Eq.  (3.35)  constitute  NN 
equations  in  NN  unknowns  for  ZK(1),  ZK(2),  ...  ZK(NN).  Upon  solving 
these  equations  and  the  ZK ( K)  determined  we  can  use  Eq.  (3.32)  to 
find  any  value  of  YM  for  a  given  value  of  UX  in  the  interval  (UX(1), 
UX(NN) ) . 

"Spline  fit  interpolation  in  the  table  of  values  (UX(1),  YM( 1 ) ) , 

( UX ( 2 ) ,  YM ( 2 ) ) ,  ....  (UX(NN),  YM(NN) )  consists  of  determining  which 


two  points  the  given  value  of  UX  lies  between  and  then  finding  YM  by 
the  formula: 

YM  *  CN(1,K)*(UX(K+1)  -  UX) 3  +  CN(2,K)*(UX  -  UX(K))3  + 

CN(3,K)*(UX(K+1)  -  UX)  +  CN(4,K)  (UX  -  UX(K) )  (3.39) 

where  the  constants  CN(1,K),  ,  CN(4,K)  can  be  determined  by  solving 

Eqs.  (3.35),  (3.37)  and  (3.38)  for  ZK(1),  ZK(2) . ZK(NN) : 

CN(1,K)  »  ZK(K)/6DK(K)  (3.40) 

CN(2,K)  =  ZK(K+1)/6DK(K)  (3.41) 

CN(3,K)  »  YM( K)/DK( K)  -  ZK(K)*DK(K)/6  (3.42) 

C N( 4 , K )  *  YM(K+1)/DK(K)  -  ZK(K+1)  *  DK(K)/6  (3.43) 


By  making  the  following  definitions 

PK(K)  =  DK(K)/6;  K  =  1,  2,  ...  NN-1 

EL(K)  =  (YM(K)  -  YM(K-l) ) / DK ( K- 1 ) ,  K  =  2 . NN 

BK ( 1 )  =  0 

BK(K)  =  EL(K+1)  -  EL ( K ) ;  K  =  2,  ...,  NN-1 
BK(NN)  =  0 


(3.44) 

(3.45) 

(3.46) 

(3.47) 

(3.48) 


Then  we  can  write  the  system  of  Eqs.  (3.35,  (3.37)  and  (3.38)  in 


matrix  form  as: 

A  Z  =  h 

where  A  is  defined  in  Fig.  3.2. 


(3.49) 


Making  use  of  the  special  nature  of  the  above  matrix  the  ZK(I) 
can  be  found  and  then  the  constant  CN(I,J)  determined  by  the  routine 
on  page  below. 

START 

INPUTUK(K) ,YM(K) ,K*1 , . . . ,NN 
DK(K)=UX(K+1 )-UX(K) 

PK(K)*DK(K)/6 

EL(K)=( YM(K+1)-YM(K) )/DK(K) 

BK(K)*EL(K)-El(K-l)  K-2,...,NN-1 
AC(1 ,2)a-l-DK( 1)/DK(2) 

AC (1,3)*  DK ( 1 ) / DK ( 2 ) 

AC(2,3)»PK(2)-PK( 1)*AC(1 ,2) 

AC(2,2)-2(PK(1)+PK(2))-PK(1)*AC(1,2) 

AC(2,3)*AC(2,3)/AC(2,2) 

BK  ( 2 )  *BK ( 2 ) /AC( 2,2) 

AC(K,K)»2(PK(K-l)+PK(K))-Pk(K-l)*AC(K-l,K) 


•  ,*•  *"•  .* 
.  *  .  '  *  ,  * 

•  *  k  •  .  -  • 


BK ( K ) =BK ( K ) -PK ( K- 1 ) *BK ( K- 1 ) 

AC(K  ,K+1) =PK( K)/AC( K,K) 

AC(NN,NN-l)=l+DK(NN-2)/DK(NN-l)+AC(NN,NN-l)*AC(NN-l,NN) 

AC( NN,NN)=-DK( NN-2)/DK( NN- 1 ) - AX ( NN,NN-1 )*AC( NN-1 ,NN) 

BK( NN)=BK(NN-2)-AC(NN,NN-l) *(BK( NN-1) 

ZK( NN)=BK( NN)/AC( NN,NN) 

ZK(K)=3K(K)-AC(K,K+1)*ZK(K+1) ;  K-M-1....2 
ZK ( 1 ) =- AC (1,2 ) *ZK (2) -AC (1,3) *ZK (3) 

C(1  ,K) ,ZK(K)/6*DK(K) 

C(2,K)=ZK(K+1)/6*DK(K) 

C(3,K)=YM(K)/DK(K)-ZK(K)*PK(K) 

C(4,K)=YM(K+1)/DK(K)-ZK(<+1)*P(K) 

STOP" 16 

In  the  MASCOR  subroutine  the  proton  mass  correction  graph  is 
interpolated  first  for  twelve  known  points  used  as  input  and  all 
interpolated  values  for  the  mass  correction  from  proton  number  22  to 
proton  number  80  as  output.  The  neutron  mass  correction  graph  Is 
interpolated  with  twelve  known  points  as  input  and  the  values  for  the 
mass  correction  for  all  neutron  numbers  from  31  to  120  as  output.  The 
output  Is  transferred  to  the  main  program  in  a  common  statement  in  the 
arrays  YC(IP)  for  neutron  number  IP,  and  YE(IP)  for  proton  number  IP. 

The  main  program  transfers  these  corrections  to  the  binding 
energy  subroutine  and  there  the  values  are  added  to  the  rest  mass 
RM(I)  of  the  nuclide  with  the  proton  and  neutron  numbers  matching  the 
index  of  the  correction  factors.  If  for  Instance  a  nuclide  (I)  has  a 
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CHAPTER  4 


DISCUSSION  AND  CONCLUSION 

The  code  developed  in  this  thesis  calculates  the  total  gamma 
energy  released  from  the  thermal  fission  of  U-235,  produces  a  prompt 
neutron  spectrum,  and  determines  the  neutron  multiplicity  for  the 
thermal  fission  of  U-235.  The  code  was  run  varying  the  adjustable 
parameter  RR  (see  Fig.  4.1)  in  order  to  determine  the  optimum  results. 
When  the  optimum  values  were  obtained,  the  spectrum  of  prompt  neutron 
energies  was  formed  (see  Fig.  4.2). 

Using  the  value  of  1.835  for  RR,  the  total  gamma  energy  is  deter 
mined  to  be  3.2285  MeV/fission,  and  the  average  energy  per  neutron  is 
calculated  to  be  1.986  MeV/neutron.  During  the  processing  of  the 
code,  it  is  interesting  to  edit  the  number  of  neutrons  released  by 
each  of  the  individual  fragments.  The  results  of  this  edit  are  found 
in  Appendix  I.  Another  advantage  of  this  code  is  that  actual  values 
of  the  energy  released  in  fission  are  computed  for  each  fragment, 
versus  the  use  of  average  values  for  a  heavy  and  light  fragment, 
thereby  increasing  the  accuracy  of  the  computed  excitation  energies. 

The  prompt  neutron  spectrum  produced  in  the  program  results  in  a 
spectrum  that  shows  an  abundance  of  very  low  energy  neutrons  in  the 
first  energy  bin  (that  is  less  than  0.5  MeV).  This  can  be  explained 
by  the  fact  that  evaporation  theory  allows  for  the  emission  of  any 
neutron  from  a  fragment  that  has  enough  residual  energy  greater  than 
the  separation  energy  plus  the  pairing  energy.  No  compensation  for 
allowed  or  not  allowed  transitions  based  on  angular  momentum 
considerations  is  accounted  for. 


Fig.  4.1  This  graph  was  used  to  determine  a  value  for  the  adjustable 
parameter  RR.  The  line  labeled  with  X's  represents  the 
neutron  multiplicity  and  units  of  the  Y-axis  are  neutrons/ 
fission.  The  line  labeled  with  Y's  represents  the  average 
neutron  energy  and  the  units  for  the  Y-axis  are  MeV/neutron 
The  line  labeled  with  Z's  is  the  total  gamma  energy  and  the 
units  for  the  Y-axis  are  MeV/flssion.  The  optimum  value  Is 
determined  to  be  RR  *  1.835. 


Fig.  4.2  This  prompt  neutron  energy  spectrum  was  generated 

from  the  code  developed  In  this  thesis.  Each  energy 
bin  is  0.5  MeV.  Each  square  represents  the  number 
of  neutrons  found  in  each  particular  energy  bin. 


The  gamma  energy  is  very  low.  As  Terrell  noted  in  his  paper,  the 
value  of  the  gamma  energy  he  obtained  is  approximately  4.9  MeV.  This 
code  obtains  a  value  of  3.2285  MeV  and  the  actual  values  based  on  a 
simplified  model  is  approximately  6.7  MeV.  This  leads  one  to  believe 
that  gamma  ray  emission  may  compete  more  strongly  with  neutron 
emission,  and  that  evaporation  theory  does  not  adequately  treat  the 
gamma  emission  problem  adequately.  Because  there  is  great  distortion 
involved  in  the  fission  process,  a  rapidly  varying  electric  field  may 
actually  induce  gamma  emission  instead  of  neutron  emission. 

Since  the  excitation  energy  is  computed  by  using  the  energy  re¬ 
lease  and  the  total  fission  fragment  kinetic  energies,  any  error  in¬ 
duced  into  these  values  can  cause  an  error  in  the  excitation  energy. 
For  example,  the  parameter  RR  is  held  constant  for  all  fission  frag¬ 
ments.  In  reality,  because  of  the  varying  sizes  of  the  fission  frag¬ 
ments,  this  parameter  will  vary,  so  the  total  kinetic  energy  will 
differ  from  the  values  compued  in  the  code. 

Another  consideration  which  can  cause  a  lower  kinetic  energy  of 
the  fission  fragments  is  that  the  nucleus  will  "cool"  with  time  so 
that  the  kinetic  energy  of  the  fragments  decrease  by  the  time  the 
second  neutron  is  emitted,  leaving  less  residual  energy  for  that  neu¬ 
tron,  and  this  can  result  in  more  residual  energy  for  gamma  emission. 

The  average  kinetic  energy  of  prompt  neutrons  emitted  in  the 
thermal  fission  of  U- 235  is  very  well  predicted  in  the  code.  Terrell 
determined  a  value  of  1.935  MeV  which  favorably  agrees  with  experi¬ 
mental  values  found  by  Cranberg  and  Leachman.  This  code  predicts  a 
value  of  1.986  MeV/neutron. 


In  summary,  this  code  was  developed  in  order  to  use  all  possible 
fission  fragment  combinations  for  the  thermal  fission  of  U-235  in 
order  to  predict  a  prompt  neutron  spectra,  average  prompt  neutrcn 
energy,  and  total  gamma  energy.  An  additional  feature  of  the  code  is 
its  ability  to  determine  the  individual  neutron  multiplicity  of  each 
fragment  produced  in  fission. 
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APPENDIX  I 


BLOCK  DATA  INPUT 


The  following  data  is  contained  in  the  code 
as  input  and  found  in  block  data  format.  It  is 
used  in  order  to  determine  the  independent  yields 
of  the  individual  nuclides  in  order  to  yield 
averge  the  neutron  multiplicity  and  total  gamma 
energy  computed  in  the  code. 
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APPENDIX  II 


COMPUTER  LISTING 


The  following  computer  listing  is  a  code  to 
determine  the  neutron  multiplicity,  total  gamma 
energy,  prompt  neutron  spectrum,  and  neutron 
kinetic  energies  for  the  thermal  fission  of 
Uranium-235. 
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APPENDIX  III 


FRAGMENT  MULTIPLICITIES 


The  following  data  is  a  compilation  of  the 
individual  neutron  multiplicities  for  each  frag¬ 
ment  produced  in  the  thermal  fission  of  Uranium- 
235.  Al.Zl  are  the  atomic  number  and  proton  num¬ 
ber  of  fragment  1,  and  A2,Z2  are  the  atomic  number 
and  proton  number  for  fragment  2.  NU  represents 
the  neutron  multiplicity  and  ihe  first  one  is  for 
fragment  1  and  the  second  one  is  for  fragment  2. 
Only  those  fragment  splits  where  one  (or  more) 
neutrons  is  emitted  by  either  or  both  fragments  is 
printed. 
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